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In this paper we explore the physics of time- dependent hydrodynamic collimation of jets 
from Young Stellar Objects (YSOs). Using parameters appropriate to YSOs we have carried 
out high resolution hydrodynamic simulations modeling the interaction of a central wind 
with an environment characterized by a toroidal density distribution which has a moderate 
opening angle of 9 P ~ 90° The results show that for all but low values of the equator to pole 
density contrast the wind/environment interaction produces strongly collimated supersonic 
jets. The jet is composed of shocked wind gas. Using analytical models of wind blown 
bubble evolution we show that the scenario studied here should be applicable to YSOs and 
can, in principle, initiate collimation on the correct scales (-R<100 AU). Comparison of our 
simulations with analytical models demonstrates that the evolution seen in the simulations 
is a mix of wind-blown bubble and jet dynamics. The simulations reveal a number of 
time-dependent non-linear features not anticipated in previous analytical studies. These 
include: a prolate wind shock; a chimney of cold swept-up ambient material dragged into 
the bubble cavity; a plug of dense material between the jet and bow shocks. We find that 
the collimation of the jet occurs through both de Laval nozzles and focusing of the wind via 
the prolate wind shock. Using an analytical model for shock focusing we demonstrate that 
a prolate wind shock can, by itself, produce highly collimated supersonic jets. 

Animations from these simulations are available over the internet at WWW address 
http: / /www. msi.umn.edu/Projects/twj/jetcol. html 



Subject headings: ISM: Jets and Outflows - hydrodynamics - star: formation 
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1. Introduction 



The propagation of jets associated with Young Stellar Objects (YSOs) has been well 
studied both analytically ( [Raga fc Kofman 1992| ) and with sophisticated numerical tools 



( Plondin, Fryxell fc Konigl 1990| , |Stone fc Norman 1994a|) . These theoretical investigations, 



in conjunction with the growing data-base of high-resolution observations, have been 
extremely useful in understanding the hydrodynamics of HH jets and HH objects. The 
origin of these jets, however, remains an issue which has yet to be resolved. The obscuring 
dust and gas surrounding young stars has made it difficult to observationally determine 
physical conditions which can constrain collimation models for either the HH jets or the 
(perhaps related) bipolar CO outflows. In the absence of these constraints a number 
of collimation mechanisms have been proposed which, broadly speaking, fall into two 
categories: pure hydrodynamic models and magnetohydrodynamic models. 

Many of the purely hydrodynamic studies produce well collimated supersonic outflows 
by invoking de Laval nozzles ( Konigl 1982 . [Raga fc Canto 1989|) . In these models an initially 



spherical stellar wind interacts with the surrounding medium and is shocked producing 
a high temperature cavity. If the walls of the cavity take on the appropriate "nozzle" 
configuration, transsonic solutions for the flow exist leading to the formation of a supersonic 
jet. There are, however, a number of problems with the de Laval nozzle scenarios. The 
nozzles in the cavity may be unstable ( |Koo fc McKee 1992[ ) and the high densities in the 



shocked gas may produce cooling distances to short to allow a "hot bubble" to form within 
the cavity flPelletier fc Pudritz 19921) . 



Another class of hydrodynamic collimation models which rely on the other extreme of 
cooling length scales was explored by Canto (1980), Canto & Rodriguez (1983) and Canto, 
Tenorio-Tagle & Rozyczka (1988). In these models strong radiative losses create a thin 
aspherical shell. After the freely expanding wind strikes the shock at an oblique angle it is 
redirected to flow along the walls of the shell. At the vertex of the aspherical (prolate) shell 
a converging conical flow is established which produces a jet. The main problem with these 
models is the size scale of the steady-state shell (based on achieving pressure equilibrium) 
which is larger than the size of observed collimation regions (R ~ 100 AU, see |Burrows~fc] 
[Stapelfeldt 19951) . 



The strong evidence supporting circumstellar disks ( |Strom 1994j ) and magnetic fields 



around T Tauri stars has led to a different set of scenarios for producing well collimated jets. 
In these MHD "disk-wind" models the outflows are centrifugally driven by a magnetized 
accretion disk. Considerable effort has gone into the theory that the magnetic field in the 
disk forces the accreting gas into co-rotation. Centrifugal acceleration and magnetic force 
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then lift the gas off the disk producing a wind which is eventually collimated into a jet as 
the field lines bend back towards the disk/star rotation axis ( [Konigl 1989| , [Pudritz 199 1| ). 
The details of the gas acceleration depend on the magnetic field configuration ( [Pelletier~& 
Pudritz 1992| , [lbdo et al 199S[ , |Wardle k Konigl 199^ , [Holland et al TM$ ), the star-disk 



interaction at the boundary layer flCamenzind 1993[ , [Najita fc Shu 19941) , and the stability 
in co-rotation ( [Uchida fc Shibata 1985|) . These disk- wind models are very promising and 
the current consensus appears to be that the outflow collimation occurs through some kind 
MHD process. However, these models also have their problems. MHD disk-wind models 
suffer from difficulties in producing the correct disk-field orientations ( [Shu 199"l~[) . Other 
magneto-gasdynamic models suffer from uncertainties in the actual field strengths and 
orientations achieved in YSOs ( |Balbus 1993| , poodman et al 1990| , |Heyer et al 1987| ). In 
addition it is not clear if or how long the field can maintain the focusing of the flow into a 
tightly collimated jet or if the collimation can be achieved on the correct scales ( [Stone fc 
Norman 19"9"415| ). 



To further complicate matters there is now increasing observational evidence that that 
collimated YSO flows are essentially unsteady. HH jets show signs of velocity variations 
( |Morse et al. 199% ) and entrainment ([Hartigan et al 1992|) , suggesting that both the driving 
of the YSO wind and its interaction with the circumstellar environment are time-dependent 
processes. 

Clearly a great deal of progress remains to be made in our understanding of the YSO 
collimation process. Because of the complexity of both the flows and the underlying physics, 
numerical simulations are an effective tool for studying outflow collimation. While there is 
an abundance of simulations of fully developed jets ( piondin, Fryxell fc Konigl 1990| , [Stone 
fe Norman 1994a]) surprisingly little numerical work has been done on their collimation 
or on the collimation of the bipolar outflows (see e.g. [Norman 1993| , Chernin & Masson 
1994). In this paper we seek to re-examine gasdynamical collimation using high- resolution 
numerical simulations. Recent numerical studies of the formation of Planetary Nebulae 
have shown that excellent collimation can be achieved through the interaction of a spherical 
central wind with a toroidal circumstellar environment ([Mellema, Eulderink fc Icke 199l| , 
Icke, Balick fc Frank 1992] , |Icke et al 1992|) . These studies have demonstrated that nonlinear 
and time- dependent gasdynamic effects provide collimating mechanisms which were 
unanticipated or not fully appreciated in previous analytical models. The time-dependent 
collimation processes have been named 'Shock- Focused Inertial Confinement' (SFIC). In the 
SFIC mechanism the interaction of the inertia of a toroidal environment with the thermal 
pressure of the shocked wind produces a well collimated jet. 



In a preliminary study [Frank fc Noriega-Crespo 1994[ (hereafter FN94) investigated 
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the SFIC mechanism in the context of YSOs. FN94 used a toroidal density environment 
which included an accretion flow. Their simulations showed that focusing at the inner shock 
could produce strongly collimated supersonic flows. While these results are promising a 
deeper understanding of the SFIC mechanism is required before a serious application to 
YSO jets can be attempted. In this paper we attempt to take some steps in this direction 
by investigating an idealized adiabatic version of SFIC collimation. Using numerical 
simulations and analytical estimates we will attempt to further identify the basic processes 
at work in SFIC collimation and put limits on their the applicability to YSOs. 

In re-examining hydrodynamic collimation our intent is not to try to push MHD models 
aside. There are many reasons for expecting magnetic fields to be important in producing 
at least some jets, particularly those associated with T Tauri stars. But a robust model of 
hydrodynamic collimation could be used to produce jets even in those cases where the MHD 
scenarios such as the disk-winds models can produce only poorly collimated winds or where 
MHD collimation is not effective on all scales relative to the full length of the jet. We note 
that the hydrodynamics we explore in this paper has elements similar to both Canto's 1980 
model and de Laval nozzles. Thus we intend to study the SFIC mechanism as a general set 
of processes which may individually operate in some form across a variety of length scales 
rather than as the definitive model for the production of YSO jets. We note here that since 
we are exploring the effect of the circum-protostellar environment on jet formation our 
models may be particularly relevent to the more deeply embedded class objects. 

The organization of the paper is as follows: In section II we describe of the numerical 
method and initial conditions used in our simulations. In Section III we provide some 
analytical estimates of the range of applicability of the SFIC mechanism with respect to 
initial conditions. In section IV we examine the results of our numerical models. In section 
V we explore the collimation mechanisms seen in the simulations. Finally in section VI we 
present our conclusions along with a discussion of some issues raised by the simulations. 



2. Numerical Method and Initial Conditions 

The numerical model (initial conditions and governing equations) we have constructed 
for our simulations capture the essential characteristics of the environment we wish to 
study: a central wind interacting with a toroidal environment. Our ultimate goal is to 
investigate the collimation of realistic YSO jets through the SFIC mechanism described in 
section I. But in this paper, we focus on the physics of SFIC collimation in an idealized 
environment. We state explicitly that the initial conditions used here are not meant to 
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be realistic in the sense of modeling an actual YSO environment. There are a number 
of scenarios for gravitational collapse that lead to star formation: collapse of a rotating 
spherical cloud ( [l'ereby, Shu fc Cassen 1954] ); collapse of a spherical cloud threaded by 
an ordered magnetic field flGalli fc Shu 199lf ); collapse of a flattened filament ( [Hartmami 



et al. 19941) . Each of the collapse scenarios listed above would produce a toroidal density 
distribution. However, the explicit form of those distributions as well as the form of the 
velocity fields they create would vary considerably from one scenario to the next. Since we 
intend to study the SFIC mechanism in the context of specific collapse scenarios in a future 
work, we use here their common characteristics as initial conditions for the environment. 

The gasdynamic interactions we wish to study are governed by the Euler equations. In 
our numerical model we express these in azimuthally symmetric cylindrical coordinates: 

<9p + ldpu L+ dpu^ = Q 
dt r dr dz 

dpu r | 1 dpu 2 | dpu T u z _ dp 
dt r dr dz dr 

dpu z + 1 dpu z u T + dpu 2 _ _dp ^ ^ 

dt r dr dz dz^ 

dE | i d(E + P )u T | g(jg+pK = (2 _ 4) 

dt r dr dz 



and 



E = \p{u r 2 + u 2 ) + J-^jP, (2-5) 



where the terms have their usual meaning. To solve these equations we use the Total 
Variation Diminishing (TVD) method of [Harten (1983)| as implemented by [Ryu et at 



(1995)| . TVD is an explicit method for solving hyperbolic systems of equations. It achieves 
second order accuracy by finding approximate solutions to the Riemann problem at each 
grid boundary while remaining non-oscillatory through the application of a lower order 
monotone scheme. The implementation of the TVD method used here is robust and 
requires even less CPU time than older methods such as the Flux Corrected Transport 
(FCT) schemes (Boris & Book 1973). 



Note that equations |2-1| through |2-5| do not include the effects of rotation and 
gravitational fields. We also leave out the effects of radiative energy losses. As said at the 
beginning of this section, we focus here on the simplest case, the purely adiabatic, purely 
gas dynamical collimation of jets. As we shall demonstrate, even under these constraints 
the flow pattern which develops is quite complex. We feel it is important to understand 
the dynamics of these flows before adding additional physics. Below we provide some 
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justification for ignoring the effects of gravity and rotation and in section 3 we discuss the 
potential role of radiative cooling in our collimation mechanisms. 

The most important feature of the environment for this study is the presence of a 
initial density contrast between pole and equator. Because gravity and rotation can both 
determine the shape of the density distribution we will use a toroidal density distribution 
with a radial power law appropriate to an in-falling cloud in the central potential of a 
protostar. The initial density distribution we chose to work with takes the following form: 

P(R,0) = " |[13P 2 (cos(0)) - 1]} (2-6) 



Note that eq |2-6| is expressed in spherical coordinates. In the rest of this paper we will use 



R and r to denote the spherical and cylindrical radii respectively where R = \/r 2 + z 2 and 
9 = tan -1 (r/z). In eq [2-6j M a is the accretion mass loss rate and M is the mass of the 
star. Equation |2-6| is a modified form of eq 96 from [Tereby, Shu fc Cassen 1984 (originally 



derived by [Ulrich 1976|) . We use it here because produces the required toroidal geometry 



as well having the R~z radial dependence, appropriate to a freely falling envelope. The 
parameter (, which determines the flattening of the cloud, is normally a function of radius 
(due to conservation of angular momentum). Since one of the principal goals of this study 
is to isolate the effect of the equator to pole density contrast (q = p e /p P ) on the SFIC 
collimation process we have modified the original equation making ( constant and treating 



it as an input parameter. Also in collapse schemes like that described [Hartmann et gl. 



1994| in which there is no rotation the pole to equator density contrast will not be a strong 
function of radius. The relation between ( and q is 

12 + 15C 



12-24C' 



(2-7) 



In the present application we set the velocity in the environment equal to zero to 
allow comparison with analytical predictions. In reality the cloud will be falling inward 



at velocities on the order of v g ~ ^JGM/R. Gravity and the infall motions of the cloud 
are, of course, the most important physical components that actually form the star. 
However, we ignore these aspects of the problem in the present study because the outward 
velocities produced by the wind/environment interaction are much greater v g and our goal 
is understanding the more restricted problem of hydrodynamic jet collimation dynamics. In 
this study we are more interested in the formation of the jets than in the formation of the 
star. Thus we can ignore the dynamical effect of the infall velocity. Although the condition 
v >> v g does not hold for all shock positions (as we will see the equatorial shock moves 
out very slowly), we found in tests that neglecting the accretion velocity does not make 
a substantial difference for the bubble structure. A similar argument can be made with 
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respect to rotational velocities since infalling material must rotate at speeds less than the 
Keplerian value. 



In Fig 1 we present a contour plot of the density distribution given by eq |2-6 . 
Superimposed on top of the density contours we plot a line marking the full-width half 
maximum of the distribution, i.e. the angle at which the density falls by half its value at 
the equator. Note that the "opening angle" of the torus defined in this way is ~ 90°. Thus 
while our initial conditions do not describe a thin disk they certainly do not correspond to 
"funnels" either and one would not, a priori, expect them to produce collimated jets. 

We assume that the temperature in the cloud is constant. The ambient pressure is then 
set by the equation of state for an ideal gas. Note that this implies an outward pressure 
gradient in the environment. But the low temperatures in the environment (T < 500 K) 
ensure small values for the sound speed and we do not see appreciable evolution of the 
environment during a simulation. 

The initial conditions for the spherically symmetric central wind are fully specified 
by its mass loss rate M w , velocity V w , and temperature T w . In the simulations the wind 
is fixed in a spherical region (R < R Q ) at the center of the grid. The use of cylindrical 
(r, z) coordinates prohibits exact specification of the inner wind "sphere" and produces a 
staircase-like pattern. This leads to small oscillations in the temperature and velocity of 
the freely flowing wind. To check that these effects do not degrade the computed solutions 
we have compared our simulations of both spherical and aspherical wind blown bubbles 
with analytical models ( [Weaver et al. 1977| , [Koo fc McKee 1992|) as well as simulations 
computed with other numerical methods in spherical coordinates ( It 1 rank 1992| ). We find 
no appreciable effects of the imperfect inner wind boundary on the bubble dynamics. The 
numerical simulations reproduced the self-similar analytical models usually to better than 
10%, but never worse than 15% in terms of shock radii and velocities. 

We have run over 30 simulations exploring a variety of initial parameters. In this paper 
we present the results of eight of these. Their input parameters are listed in Table 1. 



3. Applicability to YSOs 



The SFIC mechanism produces jets inside wind-blown bubbles with the well-known 
'three shock' structure. The first shock faces outward into the environment and is often 
called the "outer" shock. Following [Koo fc McKee 1992| we refer to this structure as the 
"ambient shock" (-R s ). The second shock faces into the central wind. It is responsible 
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for decelerating and heating the wind. In the literature it is sometimes referred to as the 
"inner" or "reverse" shock. We call it the "wind shock" (R sw ). Between these two shocks is 
the contact discontinuity (CD) separating the shocked wind gas from the shocked ambient 
gas. 

Previous simulations of SFIC collimation have shown that "flow focusing" at the wind 
shock is an important process for collimating SFIC jets. When the spherical wind strikes the 
aspherical wind shock at an oblique angle it is partly redirected into a beam aligned to the 
poles of the toroidal density distribution. Since the focusing occurs inside the wind-blown 
bubble the wind and ambient shocks must be well separated if the SFIC mechanism is 
to be effective. A similar argument holds for the presence of de Laval nozzles where the 
shocked wind fills a cavity that acts as a reservoir of thermal energy to be converted into 
bulk kinetic energy of a jet. Thus a necessary condition for the hydrodynamic collimation 
mechanisms studied here to be effective is that i? s >> R sw over at least the polar sector 
of the wind-blown bubble. Using simple analytical estimates we demonstrate below that 
it is possible to obtain this condition and, in principle, achieve hydrodynamic collimation 
through either/both de Laval nozzles and/or the SFIC mechanism on scales that are 
consistent with YSO observations. 

A shock wave is called radiative if the cooling time for the post-shock gas is shorter 
than the dynamical time scale for its evolution. If both the wind and ambient shocks are 
fully radiative then their separation will be small. Thermal energy gained at the shocks 
is quickly radiated away and, lacking pressure support, the shocked wind and shocked 
ambient material collapse on to each other forming a thin dense shell. The bubble will 
then be driven by the momentum of the wind. Given a sufficient equator to pole density 
contrast in the ambient medium, a momentum driven bubble will become highly aspherical. 
Just as in the SFIC mechanism described above, the wind in an aspherical radiative bubble 
will strike the inner shock at an oblique angle causing focusing towards the poles. In this 
situation, however, the shocked wind material must slide along the inner edge of the thin 
shell and can only form a jet directly over the poles were its conical stream converges. 
Such a configuration was the basis of Canto's 1980 model for the production of HH objects 
(see also [Tenorio-Tagle, Canto fc Rozyczka 19881 ). Numerical simulations of the SFIC 



mechanism with radiative cooling included also produce this type of flow pattern flMcllcma 
fc Frank 1996a] , |Mellema fc Frank 1996 b)) . 



When the cooling time for the shocked wind t c is comparable to, or longer than, the 
dynamical time scale for the bubble's evolution td, the post-shock wind does not lose its 
thermal energy to radiation and has enough pressure to push the wind and ambient shocks 
apart. A hot cavity of shocked wind material forms, filling a large fraction of the bubble's 
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volume. The swept-up shocked ambient material, however, remains confined to relatively 
thin shell. This is the domain where collimation through de Laval nozzles and the SFIC 
mechanism is possible. The bubble is said to be "energy conserving" or "adiabatic". 

If we wish to determine where in YSO parameter space SFIC jets might form, we must 
determine the radii at which the shocks in a wind-blown bubble begin to separate. The 
usual means of doing this is to determine the radius at which a bubble makes a transition 
from being radiative (or momentum driven) to being adiabatic (or energy driven). There 
is, however, another possibility. 



In their work on the evolution of wind blown bubbles [Koo fc McKee 1992| showed that 



between the radiative and adiabatic configurations lies another evolutionary state which 
they called the Partially Radiative Bubble (PRB). In a PRB the cooling time for the gas is 
shorter than the age of the bubble but longer than the time it takes for the unshocked wind 
to reach the wind shock i.e. t cr0 ss < ^cooi < td where t CTOSS = R s /v w . While at any time in the 
PRB stage most of the shocked wind will have cooled, the wind material which has recently 
past through the shock will still be hot enough to keep R s » R sw . Thus the appropriate 
transition radius we must find is i? t (PRB): the distance at which a radiative bubble makes 
the transition to a PRB. 

Of course not all bubbles will make a transition to the PRB stage. But as long as the 
radial density distribution is such that 

p(r) = Poir~ k , (3-1) 

the range of accretion rates M a , wind mass loss rates M w and wind velocities v w appropriate 
for YSOs, is such that all wind-blown bubbles with k = 3/2 will enter the PRB phase. 
Below we calculate the PRB transition radius _R t (PRB). 

The expansion of a radiative, momentum-driven bubble into an environment given by 
is 



eq 3-1 



-J- 2 



B.= 13 -?"•*•)" f-F, (3-2) 



127rpoi 

where _R S is the radius to the outer shock, and a = J 2/ {12 — 3k). The cooling time for post 
shock gas t c can be estimated from the familiar results of Kahn (1976) 

t c = (3-3) 

Ppre 

where C = 6 x 10~ 35 g cm -6 s 4 and p pre is the preshock density. By setting £ coo i = t cross and 
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using eqs [3-l| through |3-3| |Koo fc McKcc 1992| derived the transition time scale 



tt(PRB) 



4-k 



Substituting this into equation |3-2| gives 

i?t(PRB) 



127rp ia 2 



(3-k) 



-(21-5k)_ 



(3-4) 



AtiC 



)M w v v 



(3-5) 



Note that while t t (PRB) depends on k and poi> i?t(PRB) does not. Note also the strong 
dependence of _R t (PRB) on velocity. 

For comparison let us also calculate _R t (AD), the radius at which the bubble becomes 
fully adiabatic. [Koo fc McKee 1992| calculate the time scale for a PRB to become an 
adiabatic bubble but it depends on 7 SW the ratio of specific heats in the PRB's shocked 
wind. Since this quantity cannot be calculated in a straight forward way we will derive 
i? t (AD) by equating the age of a momentum driven bubble with the cooling time for the 
post-shock wind t c 



td- In this case eq |3-3| takes the form 

4nCviR 2 . 



t « 



(3-6) 



In eq 



we have assumed that the shock velocity v B is approximately equal to the wind 
velocity v w . Inverting eq |3-2| gives a dynamical time tdyn- Using this we find, 



Rt(AD) 



f a 



127rpoiM ft 
3- A; 



(3-7) 



In Fig 2 we show _R t (PRB) and i? t (AD) versus v w for M w = 10 -7 M Q yr _1 . The curve 
for R t (AD) was calculated using M a = 10 -6 M Q yr _1 and k = 3/2. Figure 2 demonstrates 
that the bubble enters the partially radiative stage at R < 100 AU for v w > 150 km s -1 . 
The results of |Hirth et al 1994] and recent HST images QBurrows fc Stapelfeldt 1995| ) give 
sizes for the collimation region on the order of 100 AU. In addition most HH jets are 
observed to have velocities on the order of 200 km s -1 or more. In a collimation model 
that relies on either shock focusing or de Laval nozzles these jet speeds imply even higher 
wind speeds. Thus Fig 2 shows that the shock configurations needed for hydrodynamic 
collimation are expected to begin to operate with initial parameters and on size scales 
compatible with those derived for YSOs. 
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4. 



Results 



4.1. 



Basic Flow Pattern 



In this section we focus on the results of a single jet producing simulation: case A in 
Table 1. At the end of this section we explore the role of initial conditions on the final flow. 
For now we note that while the density contrast used in case A is high (q = 70) the features 
seen are characteristic of all the other simulations where jets appear (q > 7). In Figs 3 and 
4 we present results of the case A simulation after 1035 years of evolution. Figure 3 shows 
a gray scale map of the logarithm of the density alongside of a vector map of the velocity 
field. Figure 4 shows gray scale maps of both temperature and pressure. Note that the 
darkest gray tones in Fig 3 correspond to low values of the density, whereas in Fig 4 the 
darkest gray tones correspond to high values of temperature and pressure. 

Figures 3 and 4 demonstrate that the central wind, emerging from the base of the grid, 
becomes highly focused through the interaction with the environment. While there are 
features of the overall flow pattern that resemble a wind blown bubble, the shocked wind 
has clearly been collimated into a supersonic jet. The collimation can be most clearly seen 
in the velocity vectors in Fig 3. These show a high speed flow above (behind) the wind 
shock, aligned with the z axis. 

In order to understand the nature and origin of the flow pattern we focus first on the 
density map. In the parlance of wind-blown bubble theory discussed in section 3 we can 
define the outer boundary of the interaction region or "bubble" by the shock wave driven 
into the ambient medium by the central wind. Behind this ambient shock is a shell of 
swept-up, compressed ambient gas. At any height z the highest densities in the flow are 
found in this shell. The inner boundary of the bubble is defined by the shock wave which 
faces into the central wind, decelerating and heating it. It is the mildly aspherical feature 
at the base of the computation domain surrounding the freely expanding spherical wind. 

Ignoring the flow interior to the swept-up shell for the moment, the elongation of 
the bubble can, to first order, be explained in a simple way. Note first that the ambient 
pressure can play no role in shaping the bubble. The highest pressures in environment are 
achieved in the equator where P c oc p c T c . Even there the pressure there is always orders 
of magnitude lower than the driving thermal pressure achieved in the shocked wind with 
-Psw oc p w v w 2 . Thus only the inertia of the environment affects the shape of the bubble. In 
his study of wind blown bubble dynamics Icke (1988) used Kompaneets' (1960) formalism 
to derive an expression for the evolution of the ambient shock, 



R s — R s (&, t) 



(4-1) 



13 



dt "R s 36 

where 



W + (^S 2 )} 5 (4-2) 



2 poi{d) 



Equation |4-2| shows that A = A{9) can be defined as a local acceleration parameter for the 



ambient shock. Therefore the run of A(6) determines the asphericity of the bubble. Since 
the Kompaneets approximation assumes that P sw is constant across the shocked wind cavity 
it is the environments density distribution (inertia) which determines the 9 dependence of 
A. 

The flow of shocked wind interior to the swept-up shell departs strongly from the 
expectations of wind-blown bubble theory. According to the classic theory of non-radiative 
wind-blown bubbles the flow of the hot shocked wind should be subsonic at a high uniform 
pressure P sw . But the density map in Fig 3 shows at least two sharp discontinuities at 
heights z = 8 x 10 16 cm and z = 1.35 x 10 17 cm in the shocked wind cavity. Comparison 
of the velocity, temperature and pressure maps show that these features are strong shock 
waves, which means that the flow in the cavity has been accelerated to supersonic speeds. 
Thus it is no longer appropriate to interpret the dynamics in the simulations purely in 
terms of energy driven wind-blown bubbles. Instead we have a situation which is a mix of 
jet propagation and bubble evolution physics. 

These "internal" shocks in the cavity are quite consistent with theory of supersonic 
jets. It is well known that the interaction of a jet with the surrounding medium will produce 
two shocks: a bow-shock facing into the environment; and a jet-shock facing upstream into 
the oncoming jet material QNorman 1993|) . Consideration of the pressure map in Fig 4 



demonstrates that the leading discontinuity at height z = 1.35 x 10 17 cm can be identified 
as the the jet shock. Indeed, the shock configuration is consistent with a Mach-disk as is 
expected for a terminal jet shock. The second "internal" shock wave in the body of the 
jet at z = 8 x 10 16 cm also has a Mach disk configuration. The origin of this feature is 
consistent with the crossing shocks expected in the propagation of an initially over pressured 
jet ( [Norman 19931 ). 



The bubble's ambient shock appears to double as a bow shock for the collimated 
jet. This dual-identity is another manifestation of the mix between jet and bubble 
dynamics. Examination of animations, (which can be seen at the WWW site 
http: / / www.msi.umn.edu/ Projects/twj /jetcol.html| ) , as well as plots of the various 



quantities along the axis (see Fig 5), show the region between the jet and ambient/bow 
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shocks to be more complex than might be expected for a simple jet /environment interaction. 
From the maps shown in Figs 3 and 4 one can see a structure in this region that resembles 
a kind of oblong plug not seen in previous jet simulations. We note that simulations driving 
a jet into this kind of stratified environment have yet to be performed (see |Dal Pino 1995 



for examples of jet propagation simulations in non-constant environment). We will return 
to the origin and evolution of the jet head or "plug" in the next sub-section. 

Another feature of the simulations expected from standard jet physics is the cocoon 
of "waste" material shed by the jet across the sides of Mach disk ([Norman 1993| ). This 
material is first decelerated by the jet shock and then is diverted to flow back around the 
side of the jet body. In our simulations there is an additional feature associated with the 
propagation of the jet and the cocoon. Note the presence of a relatively cool but dense 
tongues of shocked ambient material that extends from the CD at z ~ 5 x 10 16 cm into the 
shocked wind cavity. Given the cylindrical symmetry of these simulations this feature acts 
like a chimney surrounding the jet and helps to maintain its collimation. Such chimneys 
have been observed in other SFIC simulations and they appear to be an important element 
of the inertial confinement collimation mechanism ( [Mcllcma. Euldcrink fc lcke 1991 , Ickc 
Balick fc Frank 1992| , FN94). 



To assist in identifying the basic features of the simulations in Fig 5 we present cuts 
along the z axis of the density, velocity, pressure and Mach number. The wind shock 
and ambient /bow shock can be recognized in all variables at z ~ 1000 and 12500 AU 
respectively. Similarly the jet-shock and internal Mach disk are apparent at z ~ 6000 and 
9000 AU. Note that the Mach number is frame dependent and the distinction between 
subsonic and supersonic flows is sensible only in the frame of a particular shock wave. 
Detailed examination of the evolution of the wind shock shows that it progresses very 
slowly. Its rest frame and the frame defined by the stationary grid are essentially identical 
and the Mach numbers shown in Fig 5 are only correct for the flow between the wind shock 
and the internal Mach disk. 

In Fig 3 note first that the flow is clearly being accelerated from subsonic to supersonic 
velocities (M. ~ 3) after passing through the wind shock. This suggests the presence of a de 
Laval nozzle. Note also that the average density in the body of the jet, which we define to 
be the region between the wind shock and jet-shock, is < n > « 100 cm -3 which is lower 
than the density in the environment. Thus our simulations are producing light supersonic 
jets (i.e 7] = p c /pj < 1). However, since the environmental density will continue to decline 
with distance the jet will eventually become "heavy", (77 > 1), if the simulations were to be 
continued on a larger grid for a longer time. In a real protostellar environment the jet would 
probably meet the edge of the cloud before that happened and the jet would become heavy 
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more abruptly. The addition of radiative cooling will remove lateral pressure support for 
the jet and should allow it to collapse to smaller widths and higher densities ( [Raga fc Canto 



1989| ). Thus we expect the jets produced in our simulations will always become heavy at 
some point. Finally, note again the complicated structure in the "plug", the region between 
the jet-shock and the ambient/bow shock. We will return to this point in the next section. 

The continuing collimating effect of the environment can be seen by considering the 
opening angle of the jet. The opening angle of a freely expanding supersonic jet depends 
on its Mach number, 

= 2tan- 1 (-l) (4-4) 
From equation [4-4| the jet shown in Figs 3 and 4 would, if unconstrained, have an opening 



angle of at least 40°. An opening angle of <ft ~ 22° appears more appropriate to the 
simulation, indicating that continuing confinement by the chimney and the swept-up shell 
are important in the dynamics of the jet. 



4.2. Evolution 



Our simulations demonstrate that a well collimated supersonic jet develops from the 
evolution of a wind-blown bubble, the system being an interesting mix of both wind-blown 
bubble and jet dynamics. We have already identified the ambient and wind shocks 
appropriate to wind-blown bubbles and the jet shock, crossing shocks and cocoon of waste 
jet gas appropriate to jets. In order to make the evolution of these features more explicit 
we present in Fig 6 the evolution of the system through seven sequential gray scale maps of 
the density, taken every 147 years. 

There are a number of noteworthy features in Fig 6. Firstly the evolution of the wind 
shock: as the system evolves it becomes more and more aspherical, prolate geometry. 
Using the distance to wind shock in the pole R SW (P) and equator R SW (E) we can define an 
ellipticity parameter to describe its geometry, 

Detailed examination of the wind shock shows that after 200 years of evolution it assumes a 
quasi-steady configuration with an ellipticity of < e >~ .75. As we will see the asphericity 
of the wind shock plays an important role in the collimation of the flow (Sect 5.3). 
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Capturing the wind shock poses special challenges for the numerical code as it is strong 
and extends over a relatively small region. One of the disadvantages of using a cylindrical 
code is the difficulty in modeling quasi-spherical structures on a small number of grid 
points. At later times in the evolution of the models we find that numerical errors appear 
in the wind shock. These are apparent in the small "flame" shaped region of low density 
immediately behind the wind shock and close to the symmetry axis. By viewing animations 
of the simulation we have found that this feature produces a small but noticeable effect 
on the evolution of the jet. The distortion of the wind shock drives a periodic modulation 
in the post-shock velocity. The pulses can be seen in velocity plot shown in Fig 5. Fig 6 
shows that until t ~ 900 years there is a crossing shock in the jet at a distance about 
z m .5i? sw (P), which can be explained as the expansion and subsequent contraction of an 
overpressured jet ( [Norman 1993| ). The velocity pulses produced at the wind shock however 



change the crossing shock into an "internal" Mach-disk, a structure of similar character and 
origin to the internal working surfaces explored by [Biro fc Kaga 1991 . 



Fig 6 also demonstrates the role of the collimating chimney. As the system evolves 
relatively cool and dense shocked ambient gas is continually pulled off the CD. This 
material is driven upwards into the shocked wind cavity where its inertia helps maintain the 
collimation of the jet. Comparison of Figs 3 and 4 shows the correlation between pressure 
in the jet and the shape of the chimney. The kink in the chimney occurs at roughly the 
same height as the crossing shock. 

In order to test the sensitivity of the chimney structure to numerical viscosity we did a 
series of simulations with increasing resolution. We doubled the resolution from 64 x 320 
through 256 x 1280. Each grid doubling reduces the numerical viscosity by a factor of 4. 
Reducing the viscosity in this way did not effect the existence or evolution of the chimney 
other than steepening the density gradients. Due to constraints on computational time we 
have not, however, been able to continue the grid doubling and we can not at this time say 
that our simulations are fully converged. 

The physical origin of the chimney appears to be Kelvin-Helmholtz instabilities at the 
CD. Near the base of the flow a large shear gradient exists between the shocked wind and 
the shell of shocked ambient gas as can be seen in the velocity map in Fig 3. Detailed 
inspection of animations shows that the chimney develops when corrugations in the CD 
(assumed to originate from KH instabilities) are convected up by the bulk flow in the 
shocked wind cavity. 

Some aspects of the evolution of the "plug" at the head of the jet can also be followed 
in Fig 6. Consideration of this figure and density plot in Fig 5 shows that there are two 
contact discontinuities in the plug. In Fig 5 these occur at z ~ 9500 and z ~ 12000 AU 
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respectively. From Fig 5 it can be seen the contact discontinuity occurring at at z ~ 12000 
AU is the inner edge of the swept-up shell of ambient gas. The contact discontinuity at 
z ph 9500 AU marks the contact discontinuity at inner edge of what would be, in a classic 
jet, the bow shock. In our simulations however much the material between these two 
contact discontinuities originates in the stellar wind. From inspection of the early epochs of 
the simulations it appears that the plug initially forms from subsonic material injected into 
the shocked wind cavity before the jet forms. This material entered the cavity at relatively 
high densities and was then further compressed by the jet once it develops. At later times 
however Fig 3 and 6 show that additional material is added to the plug as shocked jet gas 
exiting the Mach-disk "splatters" against the shell and is driven both backward into the 
cocoon and forward into the plug. Ambient material appears to be pulled into the plug at 
these points as well during the later evolution of the jet/bubble. 

The mixture of wind-blown bubble and jet dynamics can be quantitatively explored by 
examining the evolution of three characteristic lengths: the distance to the ambient shock 
along the pole -R S (P); the distance to the ambient shock along the equator i? s (E); the radius 
of the wind shock R sw . Recall that i? sw (P) ~ _R SW (E). For a spherically expanding bubble 
both R s (t) and R sw (t) have closed form analytical expressions. 

R s (t) = X 1 [ ]rfr. (4-6) 

Poi 

flsw(t) = A 2 H 14 - (4-7) 

Poi 

where both Ai and A2 are constants of order 1. Exact expressions for these quantities can 
be found in [Koo fc McKee 1992| . We focus first on the growth of the ambient shock. In 



Fig 7 the evolution of both R S (P) and -R S (E) is plotted at 100 years intervals. In addition 
we have also plotted the growth expected for these shocks if R oc . These curves have 
been normalized to the time and distance of the first plotted point. Along the pole the 
ambient shock is clearly expanding faster than the predicted tr rate while long the equator 
it expands slower than predicted. The inner shock is also expanding more slowly than its 
predicted rate of The points plotted with an asterisk are the analytical predictions for 
the magnitudes of Rs(P), -R S (E) and R sw respectively at t — 900 years. These values were 
calculated using the appropriate mass loss rates along the equator and the pole. Recall that 
the simulations of spherical bubbles recovered both the predicted rates and magnitudes 
to within 10%. While we do not expect spherical models to recover the magnitudes of 
aspherical bubbles the growth rates should be well matched (see Pwarkadas et al 1995|) . 



From Fig 6 however it is clear that none of these quantities is recovering the analytical 
growth rates for a wind blown bubble and only R sw is within the systematic errors of the 
predicted magnitude. 
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Consideration of the ambient shock velocity along the pole V sp is also useful in 
determining the dynamics of the system. The velocity of a jet bow shock is given by the 
familiar formula 

Ks = ^et(l + V^) _1 - (4-8) 

If -R S (P) evolved solely as a jet bow shock, then 

V sp (R) = V hs (z) = V jct (l + Vv^t)- 1 (4-9) 

where 

r, = ^1. (4-10) 
Pi 

Thus the ambient shock would accelerate along the pole. If R S (P) evolved solely as wind 
blown bubble then 

V sp (R) oc [^^}\R-l (4-11) 
Poi 

and the ambient shock would decelerate along the pole. The lower right hand panel of Fig 7 
shows the actual velocity in the simulations as a function of radius. For comparison we 
have also plotted representative curves for both the evolution of a jet bow shock and a wind 
blown bubble. Apart from small variations, the velocity is roughly constant at v sp ~ 60 
km s _1 . Thus the simulations indicate that in spite of the clear presence of a jet the global 
dynamics of the system lies between that of a pressure driven bubble and supersonic jet. 



4.3. Collimation and Initial Conditions 



In order to test the limits of the hydrodynamic collimation mechanisms under study 
we have explored the parameter space of initial conditions for the simulations. We find that 
the density contrast q is the most important parameter for determining the collimation of 
the flow. In Fig 8 we show contour plots of density from four different simulations (case C, 
D, E and F) each with different initial values of the equator to pole contrast q (g = 3, 7, 14 
and 30 respectively). These plots demonstrate that collimation of the shocked wind flow 
into a jet occurs between q = 7 and q = 14. These are not extreme values. The values of q 
obtained in numerical simulations of the collapse of rotating clouds can be as high as 1000 
flYorke et al. 1993|) . Recall also that the collimation of the flow occurs without the benefit 
of additional ram pressure from the inward directed accretion flow that would occur with 
more realistic initial conditions. Since the accretion velocity goes as R^~^ the asphericity 
which develops from the density gradient will be enhanced as the equatorial shock will 
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remain at relatively smaller radii where the ram pressures are higher. The same principal 
should hold for the inclusion of gravitational potential of the protostar. 

We have also performed a set of runs to determine the effect of wind luminosity on 
the jet collimation. Using the results from [Smith et al. 19831 , |Koo fc McKee 1992| fixed 



the following limits on the effectiveness of hydrodynamic collimation in terms of the wind 
luminosity (L w = \ M W V W 2 ), 

3 < < 10 (4-12) 



In equation [4-12| L wc is a critical luminosity based on, among other things, the scale height 



of the ambient density distribution. Since no characteristic scale exists for a power-law 
distribution a direct comparison with this prediction is difficult. However, we can bracket 
the range of luminosities where the collimation seen in our simulations operates. We have 
run two simulations (case G and H) with the density contrast fixed at q = 20 but varying 
the velocity from V w = 100 km s _1 to K, = 700 km s _1 . These simulations cover a factor 



of 49 in the wind luminosity, nine times larger than the range predicted by equation [4-12 . 
We find jet collimation occurs in both simulations. Thus, given the idealizations inherent 
in our model we find no significant limits on hydrodynamic collimation over a range of 
luminosities consistent with with those observed in YSOs. In the next section we discuss 



some reasons why equation [4-12; might be wrong. 



5. Collimation Mechanisms 



From the results presented in the previous section it is clear that it is possible to 
use the interactions with the environment to produce a high degree of collimation in the 
shocked wind as well as transsonic flow. From these simulations and others studies of 
inertial collimation it appears that a number of mechanisms contribute to the production 
of supersonic jets. In this section we explore these mechanisms in more detail. 



5.1. De Laval Nozzles 



In Fig 9 we present a contour plot of model B after 300 years of evolution with cuts 
of velocity, pressure and Mach number along the pole plotted beneath it. As in Fig 5 the 
axial cuts show that immediately behind the wind shock, R SW (P), the flow passes through a 
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sonic point. We have also marked the height at which the contact discontinuity achieves a 
minimum width (ie min[rcD(<2)]) m the plots shown in Fig 8. It is clear that the constriction 
in the walls of the shocked wind cavity occurs at the same point where the flow undergoes 
a sonic transition, in other words a de Laval nozzle has formed. 

The presence of de Laval nozzles is not, in itself, surprising. As was noted in 
the introduction there is an extensive literature on de Laval nozzles as jet collimation 
mechanisms. Currently it appears that this type of mechanism is out of favor. But the 
manifestation of de Laval nozzles seen in these simulations is quite different from the 
standard steady state models found in the literature. Therefore the conclusions which led 
to the abandonment of de Laval nozzles do not apply here. 



As was discussed in the previous subsection [Smith et al. 1983 placed stringent limits on 



the parameter space of initial conditions under which stable de Laval nozzles were expected 
to form (eq [4-12] ). The theoretical basis for this conclusion was the expectation that beyond 



these limits Kelvin-Helmholtz instabilities would choke off the nozzle producing a series 
of bubbles rather than a continuous jet. The analytical arguments invoked were bolstered 
by the results of numerical simulations carried out in an earlier paper by [Norman et al. 



1981 . Comparing these papers with our results is difficult because of the different initial 
conditions and assumptions. In spite of these differences, however, a few points can be 
made. 



Firstly, and most obvious, the numerical simulations carried out by [Norman et al. 1981 



are highly underresolved by current standards (though they were state of the art at the 
time). The simulations were performed on 40 x 40 computational grids, which is almost a 
factor of 64 smaller than used in the simulations presented here. Also, one can analytically 
show that Kelvin-Helmholtz instabilities occur, but not that they will choke off the jet 
nozzle. What matters is the non-linear, time-dependent effect of the instabilities on the 
flow. We see in our simulations that the strong shear along the contact discontinuity does 
produce ripples along its surface, but rather than choke off the jet we find the instabilities 
end up aiding its collimation by providing the material for the dense chimney. 

There is another important difference between our simulations and previous studies of 
de Laval nozzles. For the most part analytical investigations have assumed a steady state 
configuration for the nozzle by matching the pressure in the bubble with the pressure in 
the ambient medium (see |Konigl 1982| for some thoughts on the evolutionary aspects of 



nozzles). In our model the bubble is over pressured and keeps expanding. Here it is the 
inertia, not the pressure of the ambient medium which provides the basis for producing a 
collimating cavity. This means that if our models are to be applied to YSOs the wind or 
the accretion or both must be taken to be time dependent allowing the configuration to 
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maintain a stable average configuration. There is considerable observational support for 
such a conclusion as the jets themselves are known to be time dependent ( [Morse et al. 
We will take up this issue again in the final section. 



Finally we note the issue of radiative time scales. De Laval nozzle models require 
the presence of high temperature gas. This has led to the conclusion that the size scales 
required, i.e. those of adiabatic bubbles, are too large to allow nozzles to collimate on the 
scales observed. However, recognition of the PRB phase discussed in Sect 3 allows de Laval 
nozzles to operate on considerably smaller scales than was previously thought (see also 
Raga fc Canto 1989| ). This implies that de Laval nozzles may still be an important aspect 



of hydrodynamic collimation mechanisms for YSO jets. 



5.2. Shock Focusing 



While de Laval Nozzles are clearly an important aspect of the collimation of jets in 
these simulations there is another mechanism at work. As was noted earlier the wind 
shock is not spherical. In almost all the simulations which produce jets we have found 
wind shock ellipticities of 0.7 < e < 0.8 (see eq |4-5|) In addition simulations of similar 
astrophysical systems (i.e. Planetary Nebulae, SN1987A, Superbubbles) have found stronger 
departures from spherical geometries with ellipticities achieved as low as e = 0.25 ([Icke| 
|1994| , |Dwarkadas et al 1995| ). When the wind shock takes on prolate geometries the radially 
streaming wind from the central source encounters it at an oblique angle. Only the normal 
component of the wind velocity will be shocked in these cases. The tangential component 
will remain unchanged. Thus the wind shock can act as a lens focusing the post-shock 
velocity vectors towards the jet axis. 

In Fig 10 we present a map of the post shock velocity vectors overlayed on contours 
of density to show the position of the inner shock and CD. In this figure the shock has 
a ellipticity of e = .79. Note that the flow vectors close to the equator emerge from 
the shock without any focusing. These gas parcels are only turned towards the axis by 
pressure gradients at some distance downstream. But as one travels up towards the pole 
the flow vectors are clearly being refracted poleward directly behind the shock. As we shall 
demonstrate below flow focusing by the wind shock can be an important component of the 
collimation process containing the potential to produce fully supersonic jets without the 
presence of a de Laval nozzle. 



-22 - 



5.3. Properties of the Inner Shock 



It is difficult to make an a priori determination of the wind shock shape based only on 
initial conditions. The degree to which the wind shock departs from spherical symmetry 
will be determined by nonlinear feedback, in terms of both thermal and ram pressures, 
from the evolving bubble. It is not yet clear how to calculate the characteristics of this 
feedback analytically. Because of this difficulty almost all analytical treatments of aspherical 
wind-blown bubble evolution have assumed the inner shock to be spherical (see e.g. |Smith| 
\et al. 1983 , |Mac-Low fc McCray 1988) . In contrast almost every study of these systems 



relying on numerical simulations has shown the inner shock to be aspherical to some degree 
( |Mac-how, McCray fc INorman 1988| , |Blondin fc Lundquist I99g ). 



Assume the shock takes on an elliptical geometry with ellipticity e defined by eq | 
To determine the degree of focusing in the post wind shock flow we must solve the oblique 
shock jump conditions. Here we repeat and extend the analysis of Icke (1988). Working in 
spherical coordinates, the angle (/5) between the ellipse and radially directed wind will be a 
function of polar angle (8) and is given by 

(3 = 6 + arctan ( ) , (5-1) 
ytan# J 

Because the wind shock is aspherical the radial distance at which the wind will encounter 
the shock depends on latitude. Thus the geometrical dilution of the wind will cause the 
pre-shock flow variables to be functions of the polar angle. We denote the pre-shock 
variables with the subscript "0" . Accounting for these variations we can use the jump 
conditions for a strong shock to express the post-shock variables (denoted with subscript 
"1") as 

vi P (9)= fc^ VoS in/3(0), (5-2) 



v lt (0)=v Q cosP(0) (5-3) 



(5-4) 



= ( 7-Hr I s (6)J-^M 2 (6)^ 2 m ~ 1 



7 + 1 



7-1 



(5-5) 
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In Eqs. 5^2 through 5^5, v\ p and V\ t are the post-shock velocity parallel and tangential 
to the shock normal, s is the sound speed and Ai is the Mach number. The total angle of 
defection (x) through the shock is given by 

2 tan ,9(0) 
(7 + 1) + (7 - 1) tan fj(9) 

In Fig 11 we show curves of x, -Pi, speed v\ = i/f i p + fit, and post-shock Mach number 
Adi as a function of polar angle for three values of the shock ellipticity e. Note that all 
points along the inner shock where x > $ have fully focused post-shock velocity vectors i.e. 
the flow at these points is directed towards the polar axis. For the three values of e shown 
all wind streamlines within 9{ = 30° of the polar axis exit the shock fully focused. Here 
we define 9f to be the latitude below which the flow is fully focused (i.e. x = m our 
simulations we found < e >~ .75 which yields 9{ = 33°. Thus, as Icke expressed in his 
original study, "even a small eccentricity causes a high degree of focusing" . 

The plots of velocity demonstrate the extent to which the wind can pass through the 
prolate inner shock and emerge into the bubble without being significantly decelerated. In 
some sense a highly prolate inner shock acts much like the reconfinement shocks explored 
by [Sanders 1983| in the context of free extragalactic jets. The plots of pressure demonstrate 



that additional collimation can be expected to occur as velocity vectors are turned poleward 
by tangential pressure gradients as is also seen in Fig 10. 

The plots of the post shock Mach number M.\ show that for more prolate shocks a 
significant region of the post-shock wind enters the bubble with supersonic velocities. Thus 
in principle it is in possible to produce a fully collimated supersonic jet purely through the 
action of shock focusing. The post shock Mach number depends upon the ratio of specific 
heats 7. Since we expect that shock focusing will begin during the Partially Radiative 
Bubble phase where the polytropic index 7 < | we have, in Fig 12, plotted max [Mi (0)] vs 7 
for 3 different ellipticities We have chosen fairly modest values of the ellipticity (e = .5, .65, 
and .8) to emphasize that supersonic flows can be expected behind a wide range of inner 
shock configurations. Fig 11 shows that as 7 decreases supersonic post shock flow can 
be achieved at relatively low ellipticity. In this case no de Laval nozzles are necessary to 
produce supersonic outflows. 



6. Conclusions 



The conclusions reached in this paper can be summarized as follows. 
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• The SFIC mechanism can, in principle, produce well collimated supersonic jets in the 
context of YSOs through purely hydrodynamic means. 

• The SFIC mechanism requires the formation of a hot shocked bubble of gas. By 
considering analytical estimates of the interaction between a fast wind and its 
surroundings it turns out that for typical YSO parameters this condition will hold 
(for outflow velocities ^ 200 km s _1 ). Although the bubble will not be fully adiabatic 
on small scales (~ 100 AU), cooling will not be efficient enough to stop the build up 
of a reservoir of hot gas. This partially radiative phase is quite important in the YSO 
case and may be in other astrophysical circumstances. A closer investigation of it 
would be interesting. 

• The flow pattern that forms in the SFIC situation is a complex mix between wind 
blown bubble physics and jet physics. While the base of the jet is a reservoir of hot, 
subsonic gas, typical of a wind-blown bubble, the jet itself is supersonic, showing 
internal shocks, and a cocoon of "waste" material. It also has the characteristic jet 
and bow shocks, in which the bow shock in fact is identical to the outer shock of the 
bubble. The jet is over pressured and less dense than its environment, but since the 
environment does not have a constant density, the jet is expected to change to a dense 
one once it has moved out to larger distances. Also in its evolution the structure 
behaves in between what is expected from analytical estimates for bubbles and jets. 

• The actual collimation in the SFIC mechanism is caused by a combination of effects. 
Firstly, the interaction with the surrounding medium creates a de Laval nozzle which 
allows a smooth transition from subsonic to supersonic flow. This nozzle is evolving 
and is a stable feature for a wide range of parameters in our simulations. We find 
none of the unstable behavior that was previously reported for these configurations. 
Secondly, the wind shock is aspherical which in itself leads to focusing of the outflow 
towards the axis. In fact, an aspherical shock may even produce super-sonic post-shock 
gas, since only the normal component of the velocity is shocked. For 7 < 5/3 this can 
happen even for mildly aspherical inner shocks. 

We again note that the models explored in this paper rely on the presence of an 
envelope of dense circum-protostellar material to produce collimated jets. Thus we feel that 
our results may be most relevent to the more deeply embedded class objects rather than 
to the more evolved objects with disks such as T Tauri stars. 

One may wonder about the long term evolution of these jets. We already pointed out 
that we expect the jets to become denser than their environment as they move out. We 
also pointed out that they are over pressured. This means that the whole bubble structure 
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expands, also laterally and that given enough time, all of the circumstellar material will be 
removed. For the region of parameter space explored here the time scale for this is ~ 10 4 
years. Given the dynamical age of some HH objects it is clear that something else must 
occur if the mechanisms discussed here are responsible for collimating the jets. However, it 
is also clear that the jets are variable in time and one can envisage a situation where the 
outflow is temporarily stopped or weakened and the ram pressure of the accreting medium 
is high enough to temporarily reverse the expansion of the bubble. The jet production 
would, therefore, also stop. This would lead to a situation similar to the one studied by 
Biro fc Raga 1994| in their numerical work on jets from time-dependent sources. Analytical 



estimates indicate that a periodic scenario of this type does indeed work ( |Mellema fc Frank] 
l996bD . 



In the SFIC model it is unavoidable that a reservoir of hot gas forms at the base of the 
jet, or even that the lower part of the jet consists of high temperature gas. This material 
would in principle emit free-free emission, observable in the radio continuum and (soft) 
X-rays. However in the deeply embedded sources we are considering these X-rays might not 
be observable. A specific comparison with observations however is premature until more 
realistic physics is added. In particular the size and physical conditions in the hot cavity 
will strongly depend on radiative cooling thus we defer specific comparisons until these 
calculations have been carried out. 

It is obvious that a substantial amount of work needs to be done before the SFIC 
mechanism can be claimed to be able to explain jets from YSOs. However given its 
efficiency and the fact that something like an interaction between in and outflows must take 
place around YSOs, we plan to explore it in some more detail in future papers. Especially 
the effects of cooling and the partially radiative bubble configuration will be studied in a 
next paper. 
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Table 1: Initial Conditions For Runs A - H 



run 


M w 




v w 


M a 




Q 


Resolution 


A 


1 x 10" 


-7 


200 


1 x 10- 


-5 


70 


256 x 1280 


B 


2 x 10" 


-7 


300 


1 x 10- 


-5 


70 


256 x 1280 


C 


1 x 10" 


-7 


200 


5 x 10" 


-6 


3 


256 x 512 


D 


1 x 10" 


-7 


200 


5 x 10" 


-6 


7 


256 x 512 


E 


1 x 10" 


-7 


200 


5 x 10- 


-6 


14 


256 x 512 


F 


1 x 10- 


-7 


200 


5 x 10- 


-6 


30 


256 x 512 


G 


1 x 10- 


-7 


100 


1 x 10- 


-5 


20 


256 x 512 


H 


1 x 10" 


-7 


700 


1 x 10" 


-5 


20 


256 x 512 
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FIGURE CAPTIONS 



Fig. 1 Initial density distribution. Shown are the log 10 contours of density from eq ^6 
with an equator to pole contrast q = 70. The two solid lines show the angle at which 
P = -5 Pmax = -5 p(90°). These occur at 6 ~ 45° making the opening angle of the 
density distribution « 90° 

Fig. 2 Transition radii as a function of velocity. Shown are the radii at which a radiative 
momentum driven bubble makes the transition to a partially radiative bubble (solid 
line) and an adiabatic bubble (dashed line). The curves are plotted for M w = 10~ 7 
M Q yr _1 , M a = 10" 6 M Q yr" 1 and M =1 M . 

Fig. 3 Density and Velocity for Model A. Shown are a gray scale map of log 10 (p) and a 
map of velocity vector field for model A after 1035 years of evolution. In the density 
map dark (light) shades correspond to low (high) densities. In the velocity field map 
vectors in the inner, freely expanding wind zone have not been plotted. Thus the first 
"shell" of vectors maps out the wind shock. 

Fig. 4 Temperature and Pressure for Model A. Shown are a gray scale map of T and 
P for model A after 1035 years of evolution. In the both maps dark (light) shades 
correspond to high (low) values. This is the reverse of the density gray scale map 
shown in Fig 3. 



Fig. 5 Axial plots for Model A. Shown are plots of log 10 (p), velocity \/v x 2 + v z 2 , pressure 
P and Mach number M. as a function of height above the equator z. All plots are 
taken after 1035 years of evolution. Each plot is an average over the first 10 zones in 
cylindrical radius r. 

Fig. 6 Evolution of density for Model A. Shown are seven gray scale maps of log 10 (p) for 
model A spaced 147 years apart. In the map dark (light) shades correspond to low 
(high) densities. Size scales are same as that shown in Figs 3 and 4. 

Fig. 7 Polar and equatorial shock evolution. In left and right upper panels the distance 
to the polar and equatorial ambient shocks obtained in the simulations is shown 
(triangles) for 9 separate times. Also shown (dashed line) are curves representing the 
growth predicted by analytical models of spherical self-similar bubbles. The point 
marked with an asterisk represents the predicted magnitude of a spherical bubble 
with identical input conditions as the simulation (along pole or equator). In the lower 
left panel an identical plot is shown for the wind shock along the equator. In the lower 
right hand panel the velocity of the polar ambient shock is presented (triangles). The 
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dashed line represents the velocity predicted for a wind blown bubble (normalized to 
the first data point). The solid line line represents the velocity evolution for a jet. 
Note that the solid line can not be properly normalized to the data but its form is 
representative of the shape of deceleration given in eq 



Fig. 8 Jet collimation and equator to pole density contrast. Shown are log 10 (p) contour 
plots of four models which differ only in the values of equator to pole density contrast 
q. Upper left: Model C q = 3. Upper right: Model D q = 7. Lower left: Model E 
q = 14. Lower right: Model F q = 3. 

Fig. 9 De Laval Nozzles and Outflow Collimation. In the upper panel a log 10 (p) contour 
map of model B after 750 years of evolution is shown. Below that are cuts are along 
of the z axis of velocity, pressure and Mach number. The points marked on each axial 
plot identify the region where the width of the channel (as measured by the contact 
discontinuity) has a minimum. This is the "throat" of the nozzle. 

Fig. 10 Shock Focusing: Model A. Shown are selected log 10 p contours identifying the wind 
shock and contact discontinuity. Also shown are velocity vectors for computational 
zones immediately downs stream of the wind shock (v < v w ). The density contours 
are log 10 p = [-20.95, -20.9, -20.85, -20.2, -20.0, -19.8, -19.6, -19.2, -19.0] 

Fig. 11 Shock Focusing: Analytical Model. Post-shock flow variables as a function of polar 
angle for 3 elliptical (prolate) shocks of differing ellipticity. These plots are for a wind 
velocity of 250 km s _1 and a wind density at the equator of 200 cm -3 . Upper left: 
total deflection angle. Upper right: Mach Number Ai. Lower left: Gas Pressure P. 
Upper right: Velocity v. The ellipticities of the shocks are e = .3 (dotted line), e = .5 
(dashed line), e = .8 (dash-dot line). In the plot of total deflection angle the solid 
line corresponds to % = 9. All points to the left of this line have post-shock velocity 
vectors that are fully focused, i.e. they point towards the z axis. 

Fig. 12 Post-Shock Mach Numbers for Non-Adiabatic Wind Shocks. Shown are the 

maximum value of the post-shock Mach number for 3 shocks of differing ellipticity as 
a function of the polytropic index 7. The ellipticities of the shocks are e = .8 (dotted 
line), e = .65 (dashed line), e = .5 (dash-dot line). 



nitia Conditions: Log_1 Density 




R (au) 



This figure "Fig3_dv.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/9606 1 42v2 



This figure "Fig4_tp.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/9606 1 42v2 



Loq 1n (Number Density) 




a. j 



Moch Number (Wind Shock Frame) 

5 F I ' ' ' ' 1 ' ' ' 




a. j 



This figure "Fig6_Tser7.jpg" is available in "jpg" format from: 



http://arXiv.org/ps/astro-ph/9606 1 42v2 



Velocity and Density 




bUU 



400 



300 



200 



. i / / / / 1 / 

t ! 

"? // t f / 

/>// / / ; 
/// / / / / 



100 







00 



200 



^ , JOO 



400 



500 



o 



1 max 



4^ 



Q 
X 



> 

5' 



CD 



! i 
i i 
i / 



O 
CO 



00 

o 

o 



I / 



I / 
I / 



Q 
O 



CD 



CJ5 



